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We first show that the currently accepted statistical mechanics for granular matter is flawed. 

The reason is that it is based on the volume function, which depends only on a minute fraction of 
all the structural degrees of freedom and is unaffected by most of the configurational microstates. 
Consequently, the commonly used partition function underestimates the entropy severely. We then 
propose a new formulation, replacing the volume function with a connectivity function that depends 
on all the structural degrees of freedom and accounts correctly for the entire entropy. We discuss the 
advantages of the new formalism and derive explicit results for two- and three-dimensional systems. 

We test the formalism by calculating the entropy of an experimental two-dimensional system, as a 
function of system size, and showing that it is an extensive variable. 

PACS numbers: 64.30.+t, 45.70.-n 45.70.Cc 


The field of granular physics is in urgent need of equa¬ 
tions of state, the traditional provider of which is statisti¬ 
cal mechanics (SM). Yet, although a granular statistical 
mechanical formalism was introduced a quarter of a cen¬ 
tury ago m, no such equations have been derived yet. 
Granular SM is entropy-based. Part of the entropy is 
structural m and corresponds to the different spatial 
arrangements of the grains, with each structural config¬ 
uration regarded as a microstate. These microstates de¬ 
pend on N s d structural degrees of freedom (DOFs) in d 
dimensions, with N s the number of contact position vec¬ 
tors (see below). The volume sub-ensemble is based on 
a volume function W, which is analogous to the Hamil¬ 
tonian in thermal SM. Namely, the probability that the 
system be at a structural microstate with volume V is 
presumed to be e~ v ' x °, in analogy to the Boltzmann 
factor e~ E ^ kBT . The factor X 0 = d(W)/dS, called the 
compactivity, is the analog of the temperature in thermal 
SM [T}{3]. Every grain configuration can support an en¬ 
semble of different boundary forces, each giving rise to a 
different internal stress microstate @H2|. The boundary 
forces, g m (to = 1, are the DOFs that determine 

the stress microstates. The combined partition function 
is 


Z = 


X0 4-0 Xy 


N a M 

14 dr n ]4 dg m 

n= 1 m—1 


(i) 


where (J.^ is the stress tensor, = Foy,- is the force 
moment tensor, and X i: j = d(Tij)/dS is the angoricity 
tensor II 0- The identity of the structural DOFs, r, 
is discussed below. The two sub-ensembles are not in¬ 
dependent [8] and the total entropy, S, is the logarithm 
of the total number of microstates, both structural and 
stress. Numerical and experimental tests of the formal¬ 
ism abound [10l415j and some inconsistencies were ob¬ 
served O’ • In particular, that the compactivity does 


not equilibrate in some systems Hz]. 

Here, we first show that this stems from a fundamen¬ 
tal problem with the formulation of the volume ensemble 
- the volume function, W in Q, is flawed in that it is 
independent of most of the structural microstates that 
it is supposed to describe. Consequently, it fails to ac¬ 
count correctly for the entire entropy. We then propose 
an improved formulation that both accounts for all the 
microstates and is amenable to analytic treatment. We 
use the new formulation to calculate the new partition 
function and the mean volume in two (2d) and three di¬ 
mensions (3d). The mean volume calculation supports 
a recent claim that an equipartition principle exists in 
these systems Humus]. The problem with the volume 
function is independent of the magnitudes of the bound¬ 
ary forces g m . Therefore, for clarity, we take these to 
be negligibly small, which allows us to neglect the force 
dependent term in Q. Including large boundary forces 
is straightforward but irrelevant for our purpose here. 

In thermal systems, the partition function is a sum over 
all microstates, each involving a unique combination of 
the values of the DOFs, giving rise to a specific value of 
the Hamiltonian, H, and hence of the Boltzmann factor. 
Therefore, H must depend on all the DOFs. If its deriva¬ 
tive with respect to any DOF vanishes identically then 
H is an incorrect measure of the system’s energy and it 
leads to miscounting of the microstates and miscalcula¬ 
tion of the entropy. Thus, dependence on all the DOFs 
is an essential test of any Hamiltonian-replacing function 
in granular SM. We demonstrate below that the volume 
function not only fails this test but it is also independent 
of almost all the structural DOFs! 

We consider an ensemble of identically-generated static 
systems in d dimensions, comprising all the mechanically 
equilibrated configurations constructed from a collection 
of iV» 1 grains. For simplicity, we constrain the mean 
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coordination number, z, to be the same for each system in 
the ensemble. Let M ~ y/~N <C N and M ~ IV 2 / 3 N 
be the number of grains in contact with the boundary 
walls in 2d and 3d, respectively. The total number of 
boundary grains is (a — 1 )M > M (a = 0(1)) and in¬ 
cludes grains that do not touch the walls (see fig. [T]). 
The connectivity is determined uniquely by the inter¬ 
granular contact position vectors and it is convenient to 
parameterize these by the vectors, r, connecting near¬ 
est contacts around grains [HI [2D]- In 2d, these vectors 
run clockwise around each grain (fig. |T]) and a similar 
parameterization exists in 3d ElUD- In both 2d and 3d, 
there are {Nz — M)/2 internal intergranular contacts and 
(Nz + M) /2 contacts altogether. 



FIG. 1. A 2d granular pack, with (a — 1 )M = 19 boundary 
grains (shaded), of which M = 10 contact the walls. The 
(solid and dashed) vectors f connect a grain’s nearest con¬ 
tacts, circulating clockwise. In 3d, the inter-contact vectors 
circulate clockwise around the facet that faces a cell, as seen 
from inside the grain. The solid vectors in the figure form a 
non-directional spanning tree: they are independent, repre¬ 
senting the independent DOFs, they reach every contact and 
the dashed vectors are linear combinations of them. There 
are aM = 29 boundary vectors, and our choice of spanning 
tree includes all of them but one. 




FIG. 2. The volume function § does not depend on the inter¬ 
contact vectors r surrounding grain A and therefore cannot 
distinguish between configurations a and b. 


To illustrate the problem with the volume function 


consider the the example in fig. 0- Its volume is 


>V = - ( | r B x f c | + I (fe + r c ) x r D | ) . (2) 


We neglect the contours of the boundary grains extend¬ 
ing outside the boundary vectors Tb-te, whose relative 
contribution to the total volume is negligible for N —> oo. 
The key point is that W does not depend on the inter¬ 
contact vectors surrounding grain A. Shifting grain A as 
in fig. Eh the volume of the system is still described by 
([2]), which depends only on the unchanged DOFs. Thus, 
dW/drA = dW/drA' = 0 and W cannot register that the 
configurations in figs. E 1 and E> are different. 

This argument is general - the volume function of any 
2d pack is (see fig. E> 


n aM -2 

w = ± y 

2 ^ 

m— 1 


2^ r k X r m+ i 
k= 1 


(3) 


where r m extends between two nearest contacts of bound¬ 
ary grain m, 1 < m < aM. This function depends only 
on the boundary contacts - it is independent of any of the 
interior ones. For N —» 00 , the boundary length scales as 
\/N while the number of configurations can be estimated 
as IV! ~ N n . In contrast, the boundary grains can make 

order V^V! ~ VN configurations. Thus, the volume 

function can register only \fN /N n = ]S[An/ 2 - n 
all the pack’s configurations - a minute fraction! This 
is equivalent to describing a gas in a room by a Hamil¬ 
tonian that depends only on the gas molecules closest 
to the walls. Clearly, such a Hamiltonian cannot ac¬ 
count for all the entropy of the system. Similarly, the 
volume function cannot be a good descriptor of the gran¬ 
ular entropy. Similarly, volume functions in 3d depend 
only on the boundary inter-contact vectors, i.e. only on 
iv( 2jv 7 / 3_Ar ) of the total number of configurations - a 
vanishingly small fraction. 

Note that this problem with the volume function is 
more basic than the recently reported failure of the uni¬ 
form measure assumption in certain systems |22j . which 
can be overcome by introducing a non-uniform measure 
in Q. Our new formulation below is similarly indepen¬ 
dent of this issue and can be used with any measure. 

Having concluded that the volume function is not a 
good equivalent of the Hamiltonian, the question is what 
could replace it. We propose a connectivity function, C, 
that does not suffer from these limitations 


Nz d 

C = 'y ^ y ' bq p . a prq a r p f} , (4) 

q,p— 1 a,/3=1 


where the sum is over all the r-vectors in the system. 
The coefficients b qp .^ a p will be identified below. The term 
W/X 0 in 0 is then replaced by C/t. We name 

r = d(C)/dS , 


(5) 
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the ‘contacture’ and it replaces the compactivity A'o- 
Here S is the entropy, i.e. the logarithm of the num¬ 
ber of all the possible configurations that the packing 
can be arranged into under the ensemble’s constraints, r 
is a measure of the connectivity fluctuations - its increase 
corresponds to more porous and less compact structures. 

To determine the coefficients b qp . at 3 , we require that C 
be additive, i.e. that the entropy of a system, made up of 
two subsystems, is the sum of their entropies. This con¬ 
strains bqp-af 3 to have no cross terms and to be indepen¬ 
dent of q and p. Additivity also constrains the connec¬ 
tivity function to be a sum over all the r-vectors, rather 
than only over an independent subset of them. We also 
require that C be independent of the coordinate system 
orientation, constraining b qp . a p to be a scalar constant 
times the unit matrix. The constant can be absorbed 
into the definition of r and we have 

Nz d 

C = Y,r q -r q = Y J R {n) -R {n) > ( 6 ) 

q—1 n—1 

where R^ = (ri Xn ,r 2Xn , ■•■) is a vector of the x n com¬ 
ponent of all the r-vectors. This connectivity-based for¬ 
mulation is not only sensitive to all the structural mi¬ 
crostates, but C also has two significant advantages: it 
has the same units in all dimensions, as the energy in 
conventional SM and unlike the volume function, and it 
is quadratic, enabling analytic calculations of the parti¬ 
tion function, as shown below. 

Expression © is not as innocuous as it may look 
since only N s — 1 of all the r-vectors are indepen¬ 
dent. We separate into three sub-vectors, R^ = 
(R-[ n \ R ^, R ^): R contains the x n component of the 
internal independent vectors and is (N s — aM )-long (see 
below); R^ contains the independent boundary contact 
vectors, of which there are aM — 1 (see below); and 
R d contains all the remaining Nd dependent vectors, 
which can be expressed in terms of R^ and R ^' 1 as: 
R^ = Ai ■ R + A 2 ■ R[ n \ where A\ and A 2 are, re¬ 
spectively, Nd x ( N s — aM) and Nd x (aM — 1) matrices. 
In terms of the independent vectors, we have 

d r 

n _ \ ' p( n ) p( n ) 1 p( n ) p( n ) 

L - 2^ Ri ■ Ri + R b ■ H b 

"=1 L (7) 

+ + A 2 R[ n) ) ■ (AyR^ + A 2 R [ n) ) . 

The independent r-vectors, of which there are many 
choices, form a (non-directional) spanning tree on the 
contact network. We constrain our choice to include the 
aM — 1 independent boundary contact vectors (fig. [l]) , 
as this makes it easier to calculate the partition func¬ 
tion. Interestingly, this number holds both in 2d and 
in 3d, which is shown as follows. In 2d, the boundary 
is a closed perimeter of aM vectors, of which aM — 1 


are clearly independent. In 3d, the boundary is a 2d 
surface, made of aM nodes and (aM/2 vectors, where 
C is the surface’s mean number of contacts per grain. 
Using Euler topological relation for planar graphs, this 
surface consists of (C/2 — 1) aM — 1 elementary loops, 
each of which has one dependent r. Thus, in 3d, there 
are QaM/2 — [(C/2 — 1) aM — 1] = aM — 1 independent 
surface vectors. 

Using the connectivity partition function is 
^ = j [R\ n) - B 1 -R[ n) +Rl n) -B 2 -Rl n) +Ri n) -B 3 -R\ n ^/r 

d 

X d N °~ aM R ( i n) d aM - 1 R { ™ ) , (8) 

n= 1 


where , B 2 , B 3 — IT A^ • A\ , 1 -f- A^ • A 2 ,2 A^ • A\ , 
respectively. The exponential makes the contribution of 
large r-vectors to the partition function negligibly small, 
allowing us to extend the integration to 00 . The contri¬ 
bution of very small r-vectors is also negligible, allowing 
us to ignore their absence. Integrating first over R .^' 1 and 
then over R^ gives the structure partition function 


Z = 


/ ( \N a -1 

[ W T ) 

[IB! II E I 


d /2 


(9) 


where E = B 2 — jB^ ■ B ^ 1 ■ B$ . The mean connectivity 
is ( C ) = T 2 d T lnZ = ( N a — l)dr/2. We see that ( C ) 
is shared amongst all the DOFs, establishing a granular 
equipartition principle similar to the one obtained in [8|. 
Explicitly, the entropy is: 


(C) 


S= ^+\nZ= - 

t 2 . 


(As—1) ln(e7TT) — ln(|Hi ||E|) (10) 


To demonstrate the use of our formalism, we analyse 
2d experimental systems, each of 1172 discs of three dif¬ 
ferent radii, produced by the 3SR Lab [53]. For each 
system we construct the contact network, choose a span¬ 
ning tree, express the dependent r-vectors in terms of 
the independent ones, calculate the matrices Ai, Bi and 
E, and then compute the entropy using eq. ( [l0| . Fig [ 3 ] 
shows the entropy for non-overlapping subsystems of dif¬ 
ferent sizes. We find that the entropy increases linearly 
with system size, i.e. it is extensive. The increase rate 
depends on the unknown r (see eq. ©). Our relation 
is linear without subtracting ln(AH), in contrast to [22] . 

Rewriting ([8]) as Z = £(1), any structural expectation 
value is (A) = C(A)/Z. E.g. the boundary vectors satisfy 




( 11 ) 


Since (rf ) is independent of system size, r is inversely 
proportional to a typical single entry of E~ l . Simi¬ 
larly, the mean magnitude square of an internal vector is 








4 



The first sum is over the boundary triangles, r n \ and r n2 
are two edges of triangle n, and p n is the vector from 
the tetrahedron apex to the contact point that these two 
edges share. The second sum is over the K n independent 
contact vectors that make p n . The angles between the tri¬ 
angle edges, a ni are distributed around 7r/3. The angles 
that the vectors r n k make with p n , cos 9 n k = \l n l\\p \ > 
are distributed around 9 = 0. Evaluating the sum by 
averaging separately over the angles and the magnitudes 
of the contact vectors, gives 


FIG. 3. Entropy vs. number of particles of the 2d experimen¬ 
tal granular systems of 1231 1. 


(' n, q -r i>q ) = rd (G 1 ) qq /2,withG = B 1 -\Bj-B 2 1 -B 3 . 
Significantly, we verified numerically that the entries of 
E _1 and G _1 are independent of system size, which es¬ 
tablishes that r is an intensive variable. 

To calculate the mean volume, we use © and define 
the interior angle between neighbour vectors q and q + 1 
along the boundary (fig. [l]) as (l - ^) 7 r + 59 q , q+1 , 
where S9 q ^ q +i is its deviation from that of a regular 
aM-sided polygon. It is straightforward to show that 
the angle between boundary vectors rj. and r m+1 is 
k (aM ~ ^ 9 , 9 + 1 )• The 2d volume is then 


2t r 


^ aM—2 m m 

= 2 Z Z r *Tm+i sin ^ — 59 q q+ i 

m—1 k— 1 q—k ' 


( 


( 12 ) 

For M 1, the sum over the constant term 2n/aM 
dominates over the fluctuations S9 qtq+ i and we take it 
out of the integral. Since k ^ m + 1, the integration over 
rkTm+i yields (|r b |) 2 = (r%) and, using (fid]), we obtain 


<Vm> 


a 2 M 2 

~r ,r 1 


S 


27r 


-U K ~ Nr 


(13) 


where Ue = Tr {i? -1 } /(aM — 1) is the average of the 
diagonal element in the matrix E~ l . Since r = 0(1) 
then (V 2 d) ~ M 2 ~ N s and the mean volume is also 
shared equally amongst the DOFs - reaffirming the gran¬ 
ular equipartition principle [8]. 

In 3d we specialise to star-like systems, where all the 
boundary contact positions are uniquely defined in terms 
of the angles from an origin in the system. The volume 
is then a sum over tetrahedra, whose apexes are at one 
of the system’s internal contacts (e.g. the closest to the 
contact network centroid) and whose bases are the trian¬ 
gular facets that make the network’s boundary 


Aft, 


V 3d = E Z 


(Tnl X r n2 ) ■ p n | = 


n—1 

-^triangles 


= E 


(r n 1 x r n2 ) 


( K n 

z 

\k =1 


(14) 


T nk 


(v 3d ) = 


-^triangles E n 
2 l T 2 


(\n |) 2 (|ri I) 


(3/2) 3 ^ 2 fVt r i an gle S dfn TT tt 1/2 3/2 

=-^- U E U G t 


(15) 


where K n is the mean number of r-vectors between the 
system centroid and the boundary triangles. Uq = 
Tr {G -1 } /(N a — aM) is the average of the diagonal el¬ 
ement in the matrix G _1 . From dimensional considera¬ 
tions, only fVtriangies ~ N 2 ^ 3 and K n ~ N 1 ^ 3 depend on 
N in (151 and hence (V 3 d) ~ N. This substantiates the 


existence of an equipartition principle in 3d too. 

To conclude, we have pointed out a flaw in the original 
entropic formalism of granular SM - the volume function 
depends only on a minute fraction of the DOFs and is in¬ 
sensitive to most microstates. This results in a significant 
underestimate of the number of microstates and hence 
of the entropy. We then proposed to replace the volume 
function by a connectivity function, which is additive and 
depends on all the structural DOFs. The compactivity 
must then be replaced by a new measure - the contacture, 
r. The new formulation was used to obtain analytical 
expressions for the entropy and for several expectation 
values, as well as to analyse 2d experimental systems. 
We verified that the entropy is extensive, r is intensive, 
and calculated the mean volume in 2d and 3d. The mean 
volume was shown to be proportional to the number of 
structural DOFs, supporting an equipartition principle 
B US USI- The flaw pointed out here probably explains 
the observations in m that the stress-based angoricity 
equilibrates in subsystems while the volume-based com¬ 
pactivity does not, casting doubt on the usefulness of 
the compactivity as a good descriptor of the structural 
fluctuations. 

It is difficult to compare our method with recent at¬ 
tempts at static granular statistical mechanics as a glass- 
like transition [24] [25] • These rely on studying the jam¬ 
ming dynamics, using conventional positions and mo¬ 
menta, energy and temperature. The lack of ergodic- 
ity forces these into questionable relations between ener¬ 
getic and structural ensembles, e.g. that each energetic 
state corresponds to one structural configuration [25] . 
The analysis of the jamming state focuses on a partic¬ 
ular, protocol-dependent state, where description based 
on force DOFs is complete H2I27 ]. At the same time 
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our SM approach addresses the full phase space, includ¬ 
ing both structural and force DOFs. 

A major advantage of the new formulation is the Gaus¬ 
sian form of the partition function in all dimensions, mak¬ 
ing possible derivation of exact results, as we demon¬ 
strated. In particular, it paves the way to an explicit 
equation of state relating the means of the volume and 
the stress. It would be interesting to revisit previous 
analyses with the new formulation, including the cou¬ 
pling between the structure and stress microstates [51128], 
and study the contacture equilibration, as in [17]. We 
look forward to numerical and experimental tests of the 
new formulation. 

This work has been funded in part by EPSRC - 
EP/H051716/1 and two Alan Howard PhD Scholarships. 
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